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Orbital-free density functional theory as an extension of traditional Thomas-Fermi theory has attracted a lot of 
interest in the past decade because of developments in both more accurate kinetic energy functionals and highly 
efficient numerical numerical methodology. In this paper, we developed a new conjugate-gradient method for the 
numerical solution of spin-dependent extended Thomas-Fermi equation by incorporating techniques previously 
used in Kohn-Sham calculations. The key ingredient of the new method is an approximate line-search scheme 
and a collective treatment of two spin densities in the case of spin-dependent ETF problem. Test calculations 
for a quartic two-dimensional quantum dot system and a three-dimensional sodium cluster,Na2i6, with a local 
pseudopotential demonstrate that the method is accurate and efficient. 

PACS numbers: 



I. INTRODUCTION 

The most attractive feature of density functional theory 
(DFT) is its use of electron density as the basic variableji 2 ^ 
which is established by the Hohenberg-Kohn theorems. 4 
The pre-modern orbital-free DFT formalism, mainly the 
Thomas-Fermi-Dirac-von Weizsacker model, 5 ' 6 - 7 ' 8 though im- 
pressive for its simplicity compared to much more com- 
plicated wave-function-based approaches, is not quite ac- 
curate. Indeed, the great success of modern DFT is due 
to the introduction of single-particle orbitals in the Kohn- 
Sham (KS) scheme, 9 which has become the mainstay of ap- 
plications of the DFT formalism. The pursuit of more ac- 
curate orbital-free density functional theory is nevertheless 
persistent, 4 . 10 ! 11 ! 12 . 13 ! 14 ! 15 . 16 ! 17 ! 18 , 19 , 20 ! 21 , 22 and has obtained new 
momentum in recent efforts in ab initio molecular dynamics 
simulations of complex systems because of its intrinsic lin- 
ear scaling behavior. z^i 2 ^ 2 ^ 2 ^ 2 ^ 9 . In this paper, We will use 
the general term, extended Thomas-Fermi (ETF), to denote all 
these approaches. 

Besides its promising power as a numerical modeling 
tool, Thomas-Fermi and its extensions also provide the start- 
ing point for a lot of other theoretical modelsi 3 ^ 3 ^ 3 ^ 3 " For 
instance, in the Strutinsky model 34 of quantum dotsp 3 ^ 
Thomas-Fermi theory accounts for classical charging effects, 
and quantum effects are incorporated by considering the 
single-particle interference in the Thomas-Fermi effective po- 
tential and the residual interaction between oscillating part of 
electron density^ 

In all these cases, an efficient numerical method to solve 
the ETF equation is a necessary ingredient. The ETF equa- 
tion can be treated either as a nonlinear self-consistent prob- 
lem, or as a constrained minimization problem. In the former 
case, a naive implementation of self-consistent iteration is not 
numerically stable because of the so-called charge-sloshing 
effect; 36 Elaborate mixing techniques have to be employed to 
achieve rapid convergence*!^ In the latter case, the simplest 
approach is the steepest-descent (SD) methodi 2 ^^ 3 ^ which, 
though simple to implement, is a well-known "poor" min- 
imization algorithm. 38 Based on analogy to the Kohn-Sham 
problem, Wang et al. formulated the energy minimization in 



terms of a damped second-order equation^ The conjugate- 
gradient (CG) method is one of the most efficient algorithms 
for numerical optimization, but its implementation to con- 
strained minimization problems is usually very complicated. 21 
Inspired by a direct minimization conjugated gradient method 
developed for the Kohn-Sham problem^ 3 ^ . in this paper, 
we develop a new CG method to solve the ETF equation, 
which is simple to implement and very efficient. 

The paper is organized as followings. In the next section, 
the ETF equation is first cast into an illuminating form, fol- 
lowed by a detailed description of the new CG method. Sec- 
tion III reports some numerical tests in a two-dimensional 
quantum dot system and a three-dimensional sodium cluster. 
The final section summarizes the paper. 



II. THEORY 

A. Formulation of the problem 

According to spin density functional theory, 12 the ground 
state total energy of a iV-electron interacting system in a lo- 
cal external potential, V cyi t{v), is a unique functional of spin 
densities, p CT (r), under the constraints, 



J p a {r)dv = N a , 



(1) 



where a — a, (3 denotes spin-up and down indices, respec- 
tively, and N a is the number of electrons of each spin, with 
N a + Np = N. The total energy can be written as (atomic 
units are used through the paper) 4 * 9 . 

Et t[p a ,Pf3] = T s [p a ,pp] + J V cxt (r)p(r)dr 

i / / p(jh^ drdr , + Exc[p ^ pph (2) 



where T s [p Qj/ o^] is the spin-dependent kinetic energy func- 
tional of the fictitious non-interacting system that has the 
same ground state spin densities as the interacting one, and 



2 



E xc [p a , pfj] is the exchange-correlation (XC) energy function- 
als. T s [p a ,pp] is related to the spin-conpensated kinetic en- 
ergy functional, T®[p], byA^L 



T s [p a ,p p ] = -T°{2p a } + -T°{2pp]. 



(3) 



In the standard Kohn-Sham scheme^^ T s ° is calculated by 
Kohn-Sham single-particle orbitals, which are eigenfunctions 
of an effective single-particle Hamiltonian. In the orbital-free 
ETF theory, however, T s ° is approximated as an explicit func- 
tional of the density 



T°{p] 



(4) 



The first term is the Thomas-Fermi kinetic energy functional, 
a local functional of density, Trp F [p] = J tTp(p(r))dr, where 
*tf (p) is the kinetic energy density of a homogeneous elec- 
tron system with electron density p, 



,5/3 



(2D) 
(3D) 



(5) 



T^[p] is the von Weizsacker functional, which is the exact 
kinetic energy functional in the limit of rapidly varying densi- 
ties. In both 2D and 3D, T^[p] takes the following forrr^ 2 ^ 



Kid - ~ / 



P(r) 



-dr. 



(6) 



A is an empirical parameter that is widely used to correct the 
over-estimation of the von Weizsacker term. 1 In this paper, 
we take A = 0.25. The third term in Eq. (0} represents 
all other non-local extensions that have been the interest of 
recent e ffortSpA2iI2iI2iI£il2il2il2i22i21 Since this paper addresses 
mainly numerical aspects, we will focus on the so-called TF- 
AW method, which neglects the third term in Eq. @. The 
techniques developed here are nevertheless general to other 
more elaborated kinetic energy functional forms. 

Instead of minimizing the total energy directly over the spin 
densities, we introduce a new quantityj2if22i22i22 ?/v(r), de- 
fined by 



p a (r) = (W(r)) 2 , 



(7) 



which can be regarded as quasi-orbitals. The total energy can 
then be written as the functional of i/j a (r), 



(8) 



and the constraint Eq. Q is transformed to the following nor- 
malization condition, 



(9) 



Taking the contraint Eq. (|9jl into account by Lagrange's mul- 
tiplier method, we define 



(10) 



a=a,0 



where fx a is the chemical potential for spin a electrons. 

The advantage of using ip a instead of spin densities as the 
basic variable is two-fold>2ii22i22i2£ (1) The requirement that 
/9o-(r) has to be positive can be cumbersome to impose in nu- 
merical implementations, which is, however, trivial in the case 
of the ipa- formulation; (2) more importantly, by introducing 
tp a we can transform the ETF problem to exactly the same 
form as the KS problem so that all efficient techniques that 
have been developed in past decades for the KS equations can 
now be easily extended to the ETF problem. The gradient of 
the total energy with respect to Vv is 



5W[ip a ,il>p] _ 5E[p ai pp] dp a {v) 



<5Vv(r) 



5p a {r) dip a (r) 



2^ CT ?/v(r) 



i 



SE[Pa,Pt3) 

(11) 



where 



K»^ cxt (r) + 



[j^-dr'+ 5E ^^ ] . (12) 
J |r-r'| $P*{r) 



For Tw with the form of Eq. we have 

6T W 



Mr) 



from which, Eq. il Q can be written concisely as 

SW 

— — - = 2X{H <T ip IT (r) - /vT/v(r)}, 
bip a (r) 

where p a = ju CT /A, and H a is defined as 



(13) 



(14) 



n^v. \P „..(>■;} ) : A-i{^^+T^(r)} (15) 



The variational principle requires 8W/8ip a (r) = 0, which 
leads to 



(16) 



which has the same form as the KS equations, but much sim- 
plified since there is only one "orbital" for each spin state. As 
in the KS problem, Eq. (I16> is a non-linear equation that has 
to be solved in a self-consistent way. 



B. A conjugate-gradient method for minimization of ETF total 

energy 

The formulation of the ETF problem using quasi-orbitals 
opens up a lot of new possibilities to solve the ETF prob- 
lem. As in the case of the KS problem, mainly there are 
two types of approaches, self-consistent and direct minimiza- 
tion methodsi^iffi Here we introduced a direct minimization 
conjugate-gradient (DMCG) method. Such a method for the 



3 



KS problem is well developed in plane wave ab initio mod- 
eling of semi-conductor material systemsp^S. It was further 
improved by Jiang et al. in their KS-DFT study of quantum 
dots3S 

Starting from an initial guess, a conjugate-gradient algo- 
rithm for a numerical optimization problem usually involves 
three steps: 36 - 38 - 4() (l) Calculate the steepest descent vector; (2) 
Construct the conjugate gradient vector; and finally (3) Up- 
date the optimization variables by moving along the conjugate 
vector direction for a certain distance that is determined either 
by an exact line search or by approximations. The complica- 
tion in the ETF problem is the normalization constraint [Eq. 
Q, or {9)]. The algorithm described below is very similar 
to what we developed for the KS problem, but there are also 
some subtle differences that are critical for optimal efficiency. 
We will describe the algorithm for the spin-dependent case, 
and its application to the spin-independent case is straightfor- 
ward. 

At the m-th iteration, the SD vector is calculated as the neg- 
ative gradient vector [Eq. ( 1141 1 (Dirac's notation for state vec- 
tors is used as in Ref |40]) 



(17) 



With ^ EE 

then calculated as 



^ = 



(m) 



^a m) ) /N a . The CG vector is 

(18) 



with 



1^ 



p(»-i> 



(19) 



for m > 1 and — for m = 1. To satisfy the normaliza- 
tion constraint of ip a , the CG vector is further orthogonalized 

(m) \ 

to tpa ) and normalized to N a , 



/(m) 



= 1- 



4 m) )(^ (ml 



'■P. 



(m) 



dm) 



ip a is then updated by 




1/2 



(20) 



(21) 



(m+l)\ _ 



(m) 



cosf 



smt 



(22) 



where the values of 9™ ln are determined by minimizing the 
total energy as a function of 9 a and 9p, 

E(9 a ,9 p ) ee E tot [p a {r;9 Q ) 7 p p {r;9^} 

ee E[i, a (v-e a )^ p {v-e p )] (23) 



with 



</v(r; 9 a ) = </£ m) (r) cos9 a + ^ n \v) sinf 



(24) 



and p a (r;8 a ) = (ip a (r', B a )) 2 , which is equivalent to min- 
imizing the Lagrangian W[ip a ,ip/3] since the normalization 
constraints are intrinsically imposed here. 

The first and second derivatives of E(8 a , Op) with respect 
to a can be obtained by 



dE{6 ai 9 p ) 
d6„ 



dr- 



SE SVv(r;0 CT ) 



'<5Vv(r; 



d9 n 



= 2\U. 



H a {B a ,9 p ) MO. 



(25) 



and 

d 2 E{9 a ,9 f3 ) 

de„de n , 



2A 



MO*) H a (9 a ,6p) MO*)) &. 



MO*) 



Me*) 



with 



H a (9 ai 9 p ) = H(r;[p a (9 a ),p p (9 p )]) (27) 
dH ^ 6 ^ -A- 1 {(24 F (2 Pff (^))W 

(28) 



d9 a , 

S 2 Exc[Pa,P/3]x . 



)MO + / dt 



r - r' 



MO*) 
MB a ) 



r) 



$P*Sp* 



99 (j 

~ 891 
_ dp a (9 a 



~M> wx9 a + ^ cos0 CT (29) 



00„ 



= 2M0*)M0*)- 



(30) 
(31) 



There are a lot of standard optimization techniques avail- 
able to solve this two-variable minimization problem (single- 
variable minimization for the spin-independent ETF). 3 ^ We 
will, however, pursue an approximate scheme similar to the 
method we used in the KS case. 40 We first transform Eq. ( I25l l 
into a more illuminating form 

dEi9 d a & f p) = 2A (-^ rinft, + 4F> ™s9 c 



xH, 



M ] cos^ + 4 m) sinfl c 



A( - A a (9 a , Bp) sm29 a + B a (9 a ,9 p ) cos2^) (32) 



with 



A 



(m) 



{9 a ,9 p )\^) 



H, 

H a {9 a , 9 j3 



(33) 



4 



and 



1>. 



(m) 



and 



where we have used the fact that both tp a 
Now we introduce the approximation 

H a (9 a ,9 fJ ) ~H„(0,0), 



so that by setting the first derivate, Eq. d32l >. to be 
obtain 



(34) 



are real. 



4,(0,0)" 



(35) 



zero, we 



(36) 



The validity of Eq. fl35i can be established by the following 
observations: The dominant parts of the effective Hamilto- 
nian H a are the kinetic energy operator and external potential, 
which are independent of 6^ [See Eg. ( 1151 : the other parts of 
H a depend on 8 a through p a {0 a ), which plays a weaker role 
in determining 6>™ m . In the Appendix, we formulate a more 
accurate approximation that goes one step beyond Ea. P5l >. in 
which the dependence of the Hatree potential on 9 a is incor- 
porated. 

The algorithm described above is very close to the band-by- 
band DMCG method used in Kohn-Sham calculationsi 36 i 39 i 40 
There is, however, one critical difference. In the spirit of 
the band-by-band scheme, one might treat two spin densi- 
ties alternatively: Iterate only one spin density at a time for 
iVband times with the other spin density fixed. Such a sequen- 
tial treatment has the disadvantage that the conjugate gradi- 
ent relaxation is disrupted every iVband times, which impairs 
partly the high-efficiency intrinsic to a conjugate-gradient al- 
gorithm. We also found in our numerical tests that the approx- 
imate line-minimization scheme using Eq. fl36l > is sometimes 
numerically unstable, and a more accurate approximation to 
0mm jjj^g tnat described in Appendix or an "exact" line search 
is required. In spite of that, the sequential scheme is still much 
more efficient than the steepest descent method; When an ex- 
act line search is required, Eq. d36l provides a very accurate 
initial guess. We will denote that treatment as the SCG (for 
sequential conjugate-gradient) method in the following sec- 
tion. In contrast, the method described above iterates the two 
spin densities simultaneously, and in doing so has taken full 
advantage of the fact that the two quasi spin orbitals are not 
required to be orthogonal to each other. In this treatment, the 
approximation made in Eq. ( I36> is found to be stable and ac- 
curate, which dramatically reduces the computation efforts; 
iterations are always carried out in the conjugate gradient di- 
rection. Because the normalization constraint is imposed in 
each iteration, there is no accumulation of numerical errors 
that may occur in some CG algorithms. We will denote this 
treatment as CCG (for concurrent conjugate-gradient). For the 
spin-independent case, these two approaches are identical. 



III. NUMERICAL TESTS 

We will report numerical results only for finite systems, 
but conclusions from these test calculations should be appli- 




FIG. 1: Comparison of approximate and exact line search meth- 
ods in spin-independent ETF. (a) 9 m i n calculated by the approximate 
method, Eq. 1361 (dot), and by exact Brent's method in a typical ETF 
calculation (line); (b) Convergence errors vs iteration number using 
approximate and exact line search methods. 



cable to periodic infinite systems as well. For other impor- 
tant components of a full implementation of the ETF theory 
to real systems, we essentially use same techniques as we 
did in the KS case: 40 We use the particle-in-the-box basis 
for the representation of quasi-orbitals, ip a , which is a vari- 
ant of plane-wave basis for finite systems; The action of the 
effective Hamiltonian on quasi-orbitals, H a t/j a is effected by 
fast-sine transform; 40 For the Hartree potential, we use the 
Fourier convolution method on a doubly extended gridi 4 ^ 4 ^ 
Local spin density approximation is used for E xc . In particu- 
lar, we use the Tanatar-Ceperley functional 46 for 2D systems 
and the Vosko-Wilk-Nusair functional 47 for 3D systems. 

We first test the performance of the new method in a 2D 
quantum dot model system with a coupled quartic oscillator 
potential (QOP)^ 

Vext(x, y)=a (^j- + by 4 - 2Xx 2 y 2 + j{x 2 y - xy 2 )r 

(37) 

with r = yjx 2 + y 2 , a = 1CT 4 , b = vr/4, A = 0.53 and 
7 = 0.2. This potential was used by the authors in numeri- 
cal studies of e-e interaction effects in quantum dotSiiSii 3 . The 
numerical results are obtained for N = 200. For the spin- 
dependent case, we consider the triplet state, i.e. N a = 101 
and N [3 = 99. 

We first test the accuracy of the approximation of Eq. 
d36i by comparing with the exact line search in the spin- 
independent ETF case. The exact 8 mul is calculated by the 
standard Brent's method, 38 in which Eq. fl36l > is used to cal- 
culate the initial guess. Fig^a) shows approximate and ex- 
act 9 min in a typical minimization process. The approximate 
0mm r e m ains in good agreement with the exact one through 
the iterations. Fig. [Qb) plots the convergence errors vs iter- 
ation number using approximate and exact line minimization 
methods, respectively. The agreement between them is almost 
perfect. 

Fig. |2 shows the convergence behaviors of three different 
minimization methods, SD, SCG and CCG, in the case of the 
spin-dependent ETF calculations. For the SD method, we use 
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FIG. 2: (Color online) Error in the total energy vs iteration num- 
ber in three different minimization methods for spin-dependent ETF: 
steepest descent (dot), sequential conjugate-gradient (dash) and col- 
lective conjugate gradient (solid). In SCG, Brent's method for line 
minimization is used, and M,and = 5. 




FIG. 3: The geometry of the Na cluster used in the test calculations 

the formalism in Ref. but with a exact line search. Us- 
ing the approximation in Eq. d36i . the computation cost for a 
single step in CCG is much less than that of SD and SCG, yet 
the convergence is still faster to attain in CCG than in SCG as 
well as in SD. The combination of these two improving fac- 
tors reduces the computation efforts in CCG by almost two 
order of magnitude with respect to SCG. 

Finally, to test the performance of our new method in more 
realistic systems, we consider a sodium cluster , Na2i6- Since 
our main purpose here is to test the performance of our new 
method, we take the geometry of the cluster as a cubic without 
relaxation, as illustrated in Fig. |3j the distance between neigh- 
boring Na atoms is taken as 4 atomic units. In this case, the 
external potential is formed by the superposition of the local 
pseudo-potential for the Na ion on different sites, 

^cxt(r) =£Vp B (|r-Ri|). (38) 

i 

For V ps (r) we use Bachelet-Hamann-Schluter's norm- 
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FIG. 4: (a) Local ion pseudopotential for Na used in test calcula- 
tions of Na2i6; (b) Error in the total energy in the Na2i6 system vs 
iteration number. 

conserving pseudopotential 49 neglecting the non-local part, 
which is illustrated in Fig. |Ha). In the test calculations, we 
use a box of 28. 3 in real space with 81 mesh points in each 
direction. We take both the change in the total energy between 
successive iterations, AE = E^ — E^ 12 ^ 1 ^, and the norm 
of the gradient as the convergence criteria. 

Fig. |3b) shows the convergence behavior of our method 
in this realistic system. As in the previous 2D model system, 
the convergence can be easily attained in about 50 iterations, 
which confirms the high efficiency of our new conjugate- 
gradient method in realistic molecular systems. 



IV. SUMMARY AND DISCUSSION 

In this paper, an efficient conjugate gradient method for di- 
rect minimization of the extended Thomas-Fermi total energy 
is proposed. The new method is inspired by a similar approach 
for the Kohn-Sham problem. The key ingredient of the new 
method is an approximate line-search scheme and a collective 
treatment of two spin densities in the case of spin-dependent 
ETF problem. The high performance of the new method was 
verified in a simple 2D model system and a realistic sodium 
cluster. 

We close the paper with two comments. First, though the 
method presented in the paper is based on transforming the 
ETF problem into a KS-like form, which is enabled by the 
presence of the von Weizsacker term, our method has more 
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general significance in terms of its mathematical structure. We 
note that the validity of the formulation in Section II. B does 
not depends on the exact form of the effective Hamiltonian, 
H a . By defining 



1 SE 

2Wa 



(39) 



our method can be easily generalized to other orbital-free DFT 
formalisms with or without the von Weizsacker term. 

Second, recently in a benchmark ETF studies on atomic 
and diatomic systems using Gaussian basis sets, Chan et 
al. found that all gradient-based methods including CG and 
quasi-Newton methods perform poorly in minimizing ETF to- 
tal energy. 22 Though we tested our new method only in the 
plane-wave type representation, the formulation of the method 
is nevertheless general, and therefore should be equally appli- 
cable to local basis represented systems. We will leave further 
investigations to future studies. 
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FIG. 5: Relati ve errors in min [Eq. IaTA using Eq. gSJ (solid) 
and Eq. <A,10> (dashed), respectively, in the spin-independent ETF 
calculation of the QOP system with N = 200. The exact 6» min is 
calculated by the standard Brent's method. 38 
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A 



(m) 



H(0) 



(m) 



0(m) 



H(0) 



= 2 



m) 


H(0) r/> 


Xi 




Ir-r'l 



A.7) 

(A.8) 

drdr', for i,j = 1,2. (A.9) 



i (m) 



APPENDIX 

The approximation in Eq. j36t can be improved by con- 
sidering the dependence of the Hartree potential on 9. To 
simplify the notation, we present the formulation only for 
the spin-independent case, and the spin index is therefore 
dropped. Instead of neglecting the dependence of H on 9 
completely [Eq.(|35j], we make a less dramatic approximation, 



H(r;9)~H(r;0) 



Ap(v';9) 



dr' 



\r — r 

where Ap(r; 9) is the change of the density 
p(r;0)-p(r;O) 
(V m) 0) cos (9 + (r) sin* 



(A.l) 



Ap(r;< 



^(™)(r)) 



with 



= — xi(r) sin 2 + %2(r) cos sin 

Xi (r) = (^^(r)) 2 - (V m )(r))' 
X2 (r) = 2^ m \v)^ m \v). 



(A.2) 

(A.3) 
(A.4) 



Using Eqs.( IA.1IA.2> . we can obtain, after some straightfor- 
ward transformation, 

A(9) = Aq-Cu sin 2 9 + C 12 sin 9 cos 9 (A.5) 
B{9) = B Q -C 12 sm 2 9 + C 2 2siii9cos9 (A.6) 



)mm j s t j len determined by requiring 



Cn — C 2 



sin A9 



C'i 



•46>}, 



) cos 29 
(A. 10) 



to be zero 
the root of Ea. dA.lQI 



Though there is no simple close expression for 
like Eq. J36I . it can be easily solved nu- 



merically using standard techniques like the Newton-Raphson 
method. 38 The increase in the computation efforts using this 
more accurate line-minimization approximation is marginal 
when compared to that using Eq.(l36>. 

We tested the new line minimization method in the QOP 
system with N = 200 in the spin-independent case. Fig. (|5jl 
shows the relative errors in 9 mm 



A r 9 r ' 



omin 
approx 



flmin 

7 cxact 



flmin 

exact 



(A.ll) 



during a typical calculation, which shows that the new line- 
minimization method based on Eq. (IA. 101 significantly im- 
proved the accuracy of the approximated 9 mm . On the other 
hand, the number of iterations required to achieve conver- 
gence remains same, which confirms further the validity of 
Ea. (l36> . We note, however, that the improvement of Eq. 
dA.10> with respect to Eq.j36> is important to in the SCG 
scheme for spin-dependent ETF calculations to guarantee the 
numerical stability, as discussed in the paper. We expect that 
the improved method will also be useful in systems where 
Eq.d3*6l might fail. 
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